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o3 ; Abstract 

£ — . ' We propose a novel method to directly determine the dark matter (DM) induced elec- 

tron/positron spectrum at the source from experimental measurements at the earth, without 
reference to specific particle physics models. The DM induced gamma rays emitted via inverse 
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Q-T Compton scattering are obtained in a model independent way. However the results depend on 

& '. 

the choice of the astrophysical electron/positron background, which is not reliably known. Nev- 
ertheless, we calculate, as an illustration, the fluxes of gamma rays from the Fornax cluster in 

>. 

the decaying DM scenario with various astrophysical backgrounds. The predictions turn out to 

(N 

be either in disagreement with or only marginally below the upper limits measured recently by 

m ' 

£ — ■ the Fermi-LAT Collaboration. This provides a strong constraint on decaying DM scenario as the 
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gamma rays may be produced in other processes besides inverse Compton scattering, such as the 
bremsstrahlung and neutral pion decays. In addition, these gamma rays in the GeV range are 



^ shown to be almost independent of choices of cosmic ray propagation model and of DM density 

profile. 
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As one of the dominant components of the universe, dark matter (DM) has yet to show its 
existence other than its gravitational effects. The nature of DM can be explored via searches 
at colliders, as well as via direct and indirect detection experiments. Recently, indirect de- 
tection of DM has attracted great attention due to the cosmic ray electron/positron excesses 

n nn 

observed by the PAMELA [1| and Fermi satellites [2|, |3[. But the interpretation of these 
experimental results is subtle. It is not easy to exclude the possibility that these excesses 
may origin from nearby astrophysical sources. Even assuming the DM annihilation/decay 
to account for the PAMELA and Fermi-LAT observations, one has to face particle physics 
model dependence. 

In this paper we will show that, it is possible to tightly constrain the decaying DM 
interpretation of electron/positron excesses in a particle physics model independent way, 
by considering the gamma rays from nearby clusters. In other words, constraints can be 
obtained without any assumptions on details of the DM model. 

Experimentally, the PAMELA and Fermi-LAT Collaborations measure only the energy 
spectra of cosmic rays at the earth. To compare with theoretical predictions, one usually 
starts from a specific DM model to calculate the fluxes at the source and their propagation 
through the Galaxy. Obviously, it is much desired to extract their fluxes at the source where 
they are generated in a model independent way. In this paper, we shall see that the e 
fluxes at the source can be obtained by solving an integral equation analytically, without 
introducing a specific DM model l . Accordingly, gamma rays emitted by these DM-induced 
energetic leptons via inverse Compton scattering (ICS) can be predicted independent of 
any DM model, which can be tested against the recent experimental measurement 5| of 
gamma rays from nearby clusters by the Fermi-LAT Collaboration. Our method should 
be applicable to both annihilating and decaying DM scenarios, but we will focus only on 
decaying DM scenario in this paper. The discussion about annihilating DM case should be 
very similar to that of decaying DM. 

Conventionally, the e propagation in the Galaxy is governed approximately by the dif- 
fusion equation 

K{E) ■ V 2 / e DM (^,r) + ± [B(E)f™(E,r)] + Q™(E,r) = . (1) 



1 Actually the same idea was discussed by Hamaguchi et al. in 2009 [J], see the "Note added" below for 
more details. 
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0.55 
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MED 
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0.0112 


4 


MAX 


0.46 


0.0765 
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TABLE I: Parameters in propagation models. MIN/MED/MAX refer to models which yield min- 
imal/medium/maximal positron flux, respectively [g|. 

Here f® M (E,r) is the DM-induced e ± number density per unit energy. K(E) stands for 
the diffusion coefficient, which can be parameterized as K(E) = K (E/GeV) a with K and 
a given in Table [I] B(E) describes the energy loss, which is effectively given as B(E) = 
E 2 /(GeV ■ te), with te = 10 16 s being a typical time scale in the Galaxy. For decaying DM 
scenario, the source term Q® M (E, r) can be expressed as 
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Here p 



DM, 



rf M , M DM and dNf M /dE are the DM density, the decay width of a particular 



decay channel, DM particle mass and the e ± spectrum per DM decay via a particular channel, 
respectively. The summation is over all possible decay channels and X(E) contains all the 
particle physics information about DM. 

Given X(E), the DM induced e^ at the earth can be determined by solving Eq. ([1]) in a 
solid flat cylinder [6|-|8| as 2 
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2 In practice, one has to truncate the infinite series to a finite sum. When E' ~ E, the series in Eq. ([3]) 
converges very slowly since there is no exponential suppression at E' = E. In this range the solution is 



better expressed in an alternative form 
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Taking the MED propagation model and Navarro-Frenk- White (NFW) DM density profile [9( as an illus- 
tration, and reordering the series in Eq. © from small to large |A mn |, we shall take the first 1413 terms 
of the series in Eq. ^ as a good approximation. This truncated sum agrees well with / e DM within 0.1% 
error in the range E' ~ E. 
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with the cylinder coordinates z G [— L, L] in the z-direction and r e [0, i?] (i? = 20 kpc) 
in radius. Here J n is the n-th order Bessel function and C„'s are successive zeros of Jo- The 
solar system is at r = 8.5 kpc. 

Surprisingly, the DM- induced e spectrum X(E) at the source can be determined in a 
DM-model independent way once f^ M (E, r Q ) is known. Eq. ([3]) is actually the so-called 
Volterra integral equation and its inverse solution can be obtained analytically as 
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K(p) = L[K(x)} 
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The source spectrum X(E) can then be determined from the DM-induced e ± at the earth 
with energies larger than E. Here L denotes the Laplace transform and L~ x its inverse. The 
L here can be performed trivially while the Cauchy's residue theorem is needed to perform 
L^ 1 analytically. This part constitutes one of the major technical hurdles of our analysis. 
We refer to the Appendix for more details about the inverse solution of Volterra integral 
equation. 

On the other hand, f® M (E,r Q ) can be obtained by subtracting off the astrophysical 
background from the e ± spectrum measured by the Fermi-LAT Collaboration. They have 



reported the e ± spectrum in the range from 7 GeV to 1 TeV 



QQ. 



We take the conventional 



"model 0" 10] of the e^ background, which can be parameterized as 11 ] 
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In practice, the infinite series will be truncated, in the same vein of Eq. ©. 
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FIG. 1: Left: fg (E,r®) extracted from the Fermi-LAT spectrum of e by subtracting off the 
background <3? ±^{E). Right: X(E) determined from f^ M (E, r©), assuming the MED propagation 
model and NFW DM density profile. 

in units of GeV _1 m _2 s _1 sr _1 with e = E/l GeV. The total e ± background flux at the Earth 
can be expressed as 



<&? ± (£, ffi ) = S [^{E) + N x ^(E) 



(9) 



with a normalization factor N. To account for the solar modulation effects, the force field 
approximation E @ = E + etfip with <f)p = 0.55 GV has been taken. In order to leave room for 
the additional DM component below 100 GeV, we choose the normalization factor N = 0.8. 
With this astrophysical e^ background, the introduction of an additional leptonic component 
from decaying DM could provide a plausible interpretation of not only Fermi-LAT e ± excess 



11 



13j). 



but also PAMELA anomaly in the positron fraction (See, e.g., 

Shown in the left part of Fig. [TJ is a fit function of f® M (E,r Q ) obtained by subtracting 
off e ± background from the Fermi-LAT data. Taking this fit function /^ )M (£ , ,r ) as an 
input, one may obtain X(E) via Eq. ([5]). Shown in the right part of Fig. [T] is the X(E) 
thus obtained for the MED propagation model and NFW DM density profile with local 
DM density p© = 0.3 GeV/cm 3 . As discussed in the Appendix, we have made certain 
approximations in obtaining X(E). To estimate the theoretical errors, we have taken X(E) 
as an input in Eq. fl3]) to get a new f® M (E,r Q ). Shown in the left part of Fig. [2] is a 
comparison of this new f® M (E, r ) with the original fit function. One sees clearly that the 
errors are very small, never beyond few percents. 

Taking this spectrum function X(E) as an input, the ICS gamma rays can be deduced 
from the scattering of energetic e ± on starlight and CMB photons. One can then check 



these predictions against experimental measurements of gamma rays from inside/outside 
the Galactic halo. We remind that the constraints obtained in this way does not depend on 
any details of the DM model. Recently, Fermi-LAT Collaboration has measured gamma rays 
from nearby clusters of galaxies with an 18-month data set [5(. These clusters are supposed 
to be highly DM dominated and isolated at high galactic latitudes. High signal-to-noise 



ratios are anticipated 
dependent studies 



or g amma-ray observations targeting nearby clusters. Recent model- 
Il5| have shown that gamma rays from the Fornax cluster provide 
the strongest constraint for decaying DM. In the following we will focus on the DM induced 
gamma rays from the Fornax cluster. Certainly, there may exist other sources in clusters 
that can emit gamma-rays, besides DM annihilation/decay Nevertheless the ICS gamma- 
rays predicted from X(E) can give theoretical lower limits on the gamma ray flux. In the 
Fornax cluster, the ICS gamma rays comes mainly from the scattering of e on the CMB 
photons, while the effects of dust and starlight can be neglected [15l.llq|. Treating the Fornax 



cluster as a point source, we follow the same method in [15] to calculate the ICS gamma rays 
semi-analytically. Shown in the right part of Fig. [2] is the predicted gamma ray spectrum, 
which seems to disagree with the Fermi-LAT measurements of gamma rays [5| in the range 
of 1 — 10 GeV. Here we have considered the uncertainties from the total DM mass of the 
Fornax cluster. The corresponding viral masses M200, M500 and their error bars are adopted 
from [171 ] . 

We now address other relevant astrophysical uncertainties. For the ICS of elec- 
tron/positron on the CMB, the astrophysical uncertainty must be small as the spectrum 
of CMB photons is well known. The main uncertainties arise from choices of propagation 
model and of DM halo profile in determination of X(E) from f® M (E,r & ). Shown in the 
left part of Fig. [3] are the X(E)'s obtained by using the MED, MIN and MAX propagation 
models respectively, with the default NFW DM density profile. Shown in the right part of 



Fig. [3] are the X(i?)'s corresponding to the NFW, Einasto 18J and Isothermal 19( density 
profiles respectively, with the MED propagation model fixed. One sees that the choice of 
DM halo profile has almost invisible impact on the determination of X(E). This is be- 
cause the energetic leptons can not propagate a long distance and different DM profiles have 
very similar behavior except for the region near the Galaxy center. The choice of propa- 
gation models do introduce large uncertainties into the determination of X(E), but only 
for energies less than about 300 GeV. This is because, very high energy leptons must come 
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FIG. 2: NFW halo profile and MED propagation model are assumed. Left: The difference between 
the red solid line and the blue dashed line can be viewed as a demonstration of the theoretical 
error in determining X(E), as explained in the text. Right: The correspondingly predicted ICS 
flux of photon is shown in the Fornax cluster. Experimental upper limits are taken from [5J. 
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FIG. 3: Astrophysical uncertainties for the determination of X(E) from /g (E,r & ). Left: NFW 
DM density profile is assumed while propagation models are varied. Right: The MED propagation 
model is assumed while DM density profiles are varied. 



from the neighborhood of the solar system. It is reasonable to expect that the propagation 
effects should not have significant uncertainties in such a small distance. Kinematically, 
the ICS gamma rays arising from the scattering on the CMB requires the initial e^ energy 



E e <; m e v/ 'E 1 /e/2 (e is the energy of CMB photon and E^ is the energy of the final ICS 
photon). This means that the final ICS gamma rays with E 1 > 1 GeV are produced from 
the initial electrons and positrons with E e > 500 GeV, which has negligible uncertainties 
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FIG. 4: Left: Different astrophysical e backgrounds as well as the Fermi-LAT data are 
shown. Right: Taking NFW halo profile and MED propagation model, X(E) is determined from 
f® M (E,r & ) for various astrophysical e ± backgrounds. 

due to the choice of propagation model. As a result, the predicted ICS gamma rays in the 
energy range of 1 — 10 GeV have very small theoretical errors. This implies that, adopting 
the conventional astrophysical background with the normalization factor N = 0.8, decaying 
DM scenario fails to account for the e ± excesses without violating gamma-ray upper limits 
of nearby clusters observed by the Fermi-LAT Collaboration. 

However our current understanding of the astrophysical e ± backgrounds is still quite lim- 
ited. Recall that the normalization factor N = 0.8 in the conventional "model 0" background 
is chosen simply to leave room for the additional DM component below 100 GeV. Actually, 
N = 1 was alread y u sed to interpret successfully low energy pre-Fermi data, such as HEAT 



in 



20 1 an d AMS-01 21]. Recently, a full Bayesian analysis based on GALPROP was presented 
22J to predict cosmic-rays self-consistently. Taking their best fit parameters, the total e ± 



spectrum is found to be harder than that of conventional "model 0" background. As shown 
in the left part of Fig. HJ the behavior of these spectra below around 100 GeV reveals po- 
tential inconsistency between the Fermi-LAT data and other e ± observations at low energy. 
As a result, the DM induced e ± spectra f® M (E,r Q ) at the earth, which can be obtained 
by subtracting off astrophysical e^ background from the Fermi-LAT data, would even turn 
negative below 100 GeV and 130 GeV for N = 1 conventional background and the best fit 
background, respectively. This feature is certainly unphysical. 

We thus focus on more energetic e^. Eq. (jSJ) tells us that the source spectrum X(E) 
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can be reconstructed from f® M (E',r G ) at the earth with E' > E. This guarantees that the 
unphysical feature of f® M (E, r Q ) at low energy would not affect the determination of X(E) 
at the higher end of the spectrum. Adopting MED propagation model and NFW DM halo 
model, we then reconstruct X(E) via Eq. ([5]) for alternative choices of e backgrounds, 
as plotted in the right part of Fig. HI As stressed before, this determination of X(E) is 
independent of any particle physics model of DM. Unsurprisingly, the inconsistency between 
the Fermi-LAT data and the N = 1 conventional e^ background (best fit e^ background) at 
low energy leads to a negative source spectrum X(E) below 220 GeV (340 GeV) during our 
reconstruction procedure. We do not show these unphysical spectra at low energy and, for 
simplicity, assume them to be vanishing in the right part of Fig. HI Fortunately, the GeV 
ICS photons are only sensitive to the initial e^ with the energy E e > 500 GeV. 

Noticed that the alternative backgrounds lead to smaller fluxes of f® M (E',r Q ) in the 
whole energy range, compared to the N = 0.8 conventional background. As a result, the 
predicted ICS fluxes of photons should become softer. Taking the Fornax cluster as a point 
source, we show in Fig [5] the predicted ICS gamma rays in the Fornax cluster. For the 
conventional "model 0" e^ background with the normalization factor N = 1, too much 
gamma-rays from the Fornax cluster are still predicted in the energy range 1-10 GeV, which 
contradicts with the Fermi-LAT point-like upper limits. For the best fit e background from 
a Bayesian analysis, decaying DM scenario survives the experimental upper limits only if the 
Fornax cluster has a small total mass Msoo- This is a strong constraint since other processes 
besides ICS, such as the bremsstrahlung of energetic e ± and 7r° decays, would also produce 
gamma rays. However these gamma rays can not be estimated in a model-independent way. 
Anyway, little room is left for these model-dependent gamma-ray fluxes from decaying DM 
in the Fornax cluster. 

In summary, we have developed a novel method to determine the DM-induced e fluxes 
at the source from the corresponding fluxes at the earth, in a DM model independent way 
by solving the Volterra integral equation. Accordingly, gamma rays emitted by these DM 
induced energetic leptons via ICS can be predicted in a model independent way. It is 
worth noticing that the DM-induced e ± fluxes at the earth are obtained by subtracting 
off the astrophysical e background from the Fermi-LAT measurements of the total flux of 
electrons and positrons. So the prediction of ICS gamma rays depends on the choice of the 
astrophysical e ± background, which is unfortunately not well determined. As an illustration, 
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FIG. 5: The predicted ICS fluxes of photons in the Fornax cluster are shown in the left (right) part 
of the figure for the N = 1 conventional e background (the best fit e background). Experimental 
upper limits are taken from 5]. 

we calculate the flux of ICS gamma rays from the Fornax cluster in the decaying DM scenario 
with different e^ backgrounds. For the conventional e ± background with the normalization 
factor N < 1, the DM-induced ICS gamma rays from the Fornax cluster are found to exceed 
the upper limits measured by the Fermi-LAT Collaboration in the energy range of 1 — 10 
GeV. Using alternatively the best fit e^ background from a Bayesian analysis, decaying 
DM scenario survives existing observations only if the Fornax cluster has a small total mass 
M500. This is still a strong constraint as the gamma rays may be produced in other processes 
besides ICS, such as the bremsstrahlung and 7r° decays. In addition, the predicted gamma 
rays with E y ^ 1 GeV are essentially independent of choices of propagation model and of 
DM density profile. 

Note added: Two months after the first version of this paper appeared on arXiv, we 
noticed accidentally that almost the same method was proposed already in j4| to reconstruct 
the electron/positron source spectrum from the experimental fluxes at the earth. Neverthe- 
less, we used this method to calculate the ICS gamma rays from the Fornax cluster, which 
shows possible contradiction with or strong constraint from the Fermi-LAT measurements. 
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Appendix A: The Volterra integral Equation 

The Volterra integral equation is given as 

K(x-t)y{t)dt = f{x), (Al) 

with boundary condition K(0) = 1 and f(a) = 0. The solution of y(x) can be represented 

as 

, , df(x) r , „, .df(t) 

y(x) = -^ + dt R(x - t) J 
ax J a 

with 



R(x) = L l 



dx 
1 



dt 



(A2) 



- 1 



and K(p) = L [K(x)\ . 



(A3) 



pK(p) 

Here L denotes the Laplace transform and L^ 1 its inverse. As a deformation of the Volterra 
integral equation, Eq. (J3J) can be easily obtained from Eq. (1A2|) by a suitable change of 
variables. 

The kernel function in our case is K(x) 
ingly 
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where the Cauchy's Residue Theorem has been applied in the second line of the above 
equation, with Res denoting the residue. The summation Yl is over an singularities pi in 



the left half complex plane. Defining ipijp) = ^2 



Bn 
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P Amn 



the singularities in Eq. 



correspond to the zeros of ip(p) except p = 0. In practice, the infinite summation in ip{p) 
should be truncated, in the same vein of Eq. ([3]). Then the number of zeros of ip(p) 
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equals to the number of terms in the truncated series (e.g., as large as 1413 for the case of 
MED propagation model and NFW DM density profile). We have checked numerically that 
the truncated terms with smaller A mn s have negligible effects on the position of relevant 
zeros. In principle, one can find all singularities and their residues in Eq. (1A4J) . But this 
demands excessive amount of computer power and it is unnecessary, as we will see presently. 
Notice that all singularities have negative real parts as X mn < 0. Furthermore, there is 
an exponential suppression factor exp(pix) in the residue, which naturally offers a good 
damping factor. Therefore singularities in the region —200 < Re(p) < should yield a very 
good approximation, as evidenced by the left part of Fig. Notice also that the terms being 
truncated in ip{p) contribute no additional singularity in the region —200 < Re(p) < 0, as 
guaranteed by Rouche's theorem in complex analysis. So the errors of our method are well 
controlled. 

To find the roots of ip(p) = quickly, we apply the argument principle in complex analysis: 

±Jf^ dp = V-P, (A5) 

2m J c ip(p) 

if ip{p) is a meromorphic function inside and on some closed contour C and have no zeros 

or poles on C N and P denotes the number of zeros and poles of ip(p) inside the contour 

C respectively. Each zero is counted as many times as its multiplicity while each pole is 

counted as many times as its order. So we divide the region —200 < Re(p) < into 

many small regions and do such a integral around the contour of each small region to learn 

the distribution of the zeros of ip(p). Finally we use Newton's method to locate the zeros 

accurately in the corresponding small regions. This method is especially useful to locate the 

zeros off the real axis. 
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